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continuous valued, we develop quantum particle filters as a practical computational method for 
quantum parameter estimation. 
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I. INTRODUCTION 

Determining unknown values of parameters from noisy 
measurements is a ubiquitous problem in physics and en- 
gineering. In quantum mechanics, the single-parameter 
problem is posed as determining a coupling param- 
eter ^ that controls the evolution of a probe quan- 
tum system via a Hamiltonian of the form i?^ = ^i/o 
[ll|2j[3l|4l[5l|6j[7]. Traditionally, an estimation procedure 
proceeds by (i) preparing an ensemble of probe systems, 
either independently or jointly; (ii) evolving the ensem- 
ble under H^; (iii) measuring an appropriate observable 
in order to infer ^. The quantum Cramer-Rao bound 
[TJ [21 [51 [SJ [H] gives the optimal sensitivity for any possi- 
ble estimator and much research has focused on achieving 
this bound in practice, using entangled probe states and 
nonlinear probe Hamiltonians [HI 113 EI] • 

Yet, it is often technically difficult to prepare the exotic 
states and Hamiltonians needed for improved sensitivity. 
Instead, an experiment is usually repeated many times 
to build up sufficient statistics for the estimator. In con- 
trast, the burgeoning field of continuous quantum mea- 
surement |12j provides an opportunity for on-line single- 
shot parameter estimation, in which an estimate is pro- 
vided in near real-time using a measurement trajectory 
from a single probe system. Parameter estimation via 
continuous measurement has been previously studied in 
the context of force estimation [13J and magnetometry 
[T¥j . Although Verstraete et. al develop a general frame- 
work for quantum parameter estimation, both of [T31 [T?] 
focus on the readily tractable case when the dynamical 
equations are linear and the quantum states have Gaus- 
sian statistics. In this case, the optimal estimator is the 
quantum analog of the classical Kalman filter P^SIITCI [Tfj. 

In this paper, we develop on-line estimators for contin- 
uous measurement when the dynamics and states are not 
restricted. Rather than focusing on fundamental quan- 
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tum limits, we instead consider the more basic problem 
of developing an actual parameter filter for use with con- 
tinuous quantum measurements. By embedding parame- 
ter estimation in the standard quantum filtering formal- 
ism 12], we construct the optimal Bayesian estimator 
for parameters drawn from a finite dimensional set. The 
resulting filter is a generalized form of one derived by 
Jacobs for binary state discrimination [18]. Using re- 
cent stability results of van Handel [T^], we give a simple 
check for whether the estimator can successfully track to 
the true parameter value in an asymptotic time limit. 
For cases when the parameter is continuous valued, we 
develop quantum particle filters as a practical computa- 
tional method for quantum parameter estimation. These 
are analogous to, and inspired by, particle filtering meth- 
ods that have had much success in classical filtering the- 
ory [20, ,21^ . Although the quantum particle filter is nec- 
essarily sub-optimal, we present numerical simulations 
which suggest they perform well in practice. Through- 
out, we demonstrate our techniques using a single qubit 
magnetometer . 

The remainder of the paper is organized as follows. 
Section |ll] reviews quantum filtering theory. Section III 
develops the estimator and stability results for a param- 
eter from a finite-dimensional set. Section |IV| presents 
the quantum particle filtering algorithm, which is appro- 
priate for estimation of a continuous valued parameters. 
Section W\ concludes. 



II. QUANTUM FILTERING 

In this section, we review the notation and features of 
quantum filtering and quantum stochastic calculus, pre- 
dominantly summarizing the presentation in |12j . which 
provides a more complete introduction. In the general 
quantum filtering problem, we consider a continuous- 
stream of probe quantum systems interacting with a tar- 
get quantum system. The probes are subsequently mea- 
sured and provide a continuous stream of measurement 
outcomes. The task of quantum filtering is to provide an 
estimate of the state of the target system given these in- 
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direct measurements. In the quantum optics setting, the 
target system is usually a collection of atomic systems, 
with Hilbert space Ha and associated space of opera- 
tors A. The probe is taken to be a single mode of the 
quantum electromagnetic field, from which vacuum fluc- 
tuations give rise to white noise statistics. 

In the limit of weak atom-field coupling, the joint 
atom-field evolution is described by the following quan- 
tum stochastic differential equation (QSDE) 

dUt = (^LdAl - LUAt - ^L^Ldt - iHd?j Ut, (1) 

where L g -4 is an atomic operator that describes the 
atom-field interaction and H € A is the atomic Hamil- 
tonian. The interaction-picture field operators dAt^dA^ 
are quantum white noise processes with a single non-zero 
Ito product dAtdAl — dt. 

For any atomic observable Xa G A, the Heisen- 
berg evolution or quantum flow is deflned as (^a)* = 
jt{XA) = U}{Xa I)Ut- Apphcation of the Ito rules 
gives the time evolution as 

djtiXA) - jtiC[XA])dt+jti[L^,XA])dAt+jt{[XA,L])dAj 

(2) 

with Lindblad generator 
C[XA] = i[H,XA]+L^XAL-^L^LXA-^XAL^L. (3) 

Similarly, the observation process, which we take to be 
homodyne detection of the scattered field, is given by 
Mt = U}{At + Al)Ut- The Ito rules give the correspond- 
ing time evolution 

dMt = jt{L + L^)dt + dAt + dA\. (4) 

Together, Q and Q are the system-observation pair 
which define the filtering problem. The quantum flow de- 
scribes our knowledge of how atomic observables evolve 
exactly under the joint propagator in (jlj, but it is 
inaccessible since the system is not directly observed. 
Nonetheless, the scattered fields as measured in Q carry 
information about the atomic system, providing a contin- 
uous measurement of the observable L -\- L\ albeit cor- 
rupted by quantum noise. The quantum filtering problem 
is to find tt^-'^a] = IE(jt(XA)|M[o (j), the best estimate 
(in a least squares sense) of an atomic observable condi- 
tioned on the measurement record. We invite the reader 
to consult [m 122] for details on deriving the recursive 
form of this filter, which is governed by the (classical) 
stochastic differential equation (SDE) 

dTTt[XA] = Tit[C[XA]]dt + {-KtlL^XA + XaL] 

-T:t[L^ + L\T:t[XA]) x [dMt - T:t[L + L^]dt). (5) 

We see that this is an entirely classical filter, driven by 
the classical measurement stream Yj. Oftentimes, it is 
more convenient to work with the adjoint form of the 



equation for pt, which satisfies Tr[XAPt] = 7ri[XA] for 
all Xa G a. The state pt is often called the conditional 
density matrix. The SDE or stochastic master equation 
(SME) for Pt is then 

dpt = -i[H,pt]dt+{LptL^ - ~ \ptL^L)dt 

+ (Lpt + ptL^ - Tr[(i + L^)pt]pt)dWt (6) 

where the innovations process, dWt = dMt - Tr[(L + 
L^)p]dt, is a Wiener process that satisfies E[(iVF(] ~ 
and Ito rule (dWt)^ = dt. 

Quhit Example 




FIG. 1: Continuous-measurement of single qubit precessing 
in an external magnetic field 

Consider the setup depicted in Figure [T] A qubit, ini- 
tially in the pure state |-l-a;), precesses about a magnetic 
field B while undergoing a continuous measurement along 
z. In terms of the general framework, H = Bay and 
L = ^/Raz, where ^JH is the continuous measurement 
strength in the weak coupling limit. We will not dwell on 
the underlying physical mechanism which gives rise to the 
az measurement, though continuous polarimetry mea- 
surements could suffice [23^ . Plugging into ([5| , the quan- 
tum filter for the Bloch vector nt — {■nt['^x\, T^t['^y\,''^t['^z\) 
is 

dTTt[(Tx] = 2Bnt[az]dt - 2niTt[ax\dt - 2^/KTTt[ax]TTt[az\dWt 

(7) 

d-Kt[ay\ = ~2MTTt[ay\dt - 2^'Kt[ax\TTt[ay\dWt (8) 
dnt[az\ = -2BTTt[a^]dt + 2^{1 - 'Kt[azY)dWt (9) 

with innovations dWt = dMt — 2^/Knt[az]dt. It is not 
difficult to verify that the quantum filter maintains pure 
states and that the initial state uq — (1, 0, 0) remains on 
the Bloch circle in the x-z plane. Letting 6 be the angle 
from the positive a;-axis such that tan^ — TTt[az]/T^t[o'x]i 
we then simplify the filter to 

ddt = -2Bdt + Ksm{2et)dt + 2y/^cos{9t)dWt (10) 

where now dWt — dMt — 2^/11 sin Bt- Figure |2] shows a 
computer simulation of a typical measurement trajectory 
and filtered Bloch vector values when B = 0. Note that 
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FIG. 2: (Bottom) Simulated typical measurement trajectory 
for continuous Z measurement, k — 1, B = (Top) Filtered 
values of 7rt[cr:c] and TTtluz] for simulated trajectory 



in the absence of a magnetic field, the steady-states are 
the eigenstates, which are reached with probabihties 
given by the Born rule |24j . 



III. ESTIMATION OF A PARAMETER FROM A 
FINITE SET 

Using the atomic system as a probe for the unknown 
parameter ^ [29], we set the atomic Hamiltonian of the 
quantum filter to 



H = ^Ho, Ho € A. 



(11) 



Supposing we knew the true value of the parameter, the 
quantum filtering equations would give us the best least- 
squares estimate of the atomic system conditioned on the 
measurements and the knowledge of dynamics induced by 
£, through H. But given the optimality of the filter, we 
could equally well embed the parameter ^ as a diagonal 
operator S acting on an auxiliary quantum space, after 
which the filter still gives the best estimate of both system 
and auxihary space operators. Finding the best estimate 
of ^ conditioned on the measurements simply corresponds 
to integrating the equations for tt^ [S] . 

More precisely, extend the atomic Hilbert space Ha 
H^(^Ha and the operator space A i— > 'Z){H^)(^A, where 
13(7^^) is the set of diagonal operators on Ti^. Assuming ^ 
takes on N possible values {^i, . . . , £,n}, dimS(7i^) = N. 
Introduce the diagonal operator 



N 



(12) 



so that S|^i) = with G Ti.^. This allows one to 



generalize ( 11 ) as 



H^E(g)HoeTi{n^)(g>A. 



(13) 



Any remaining atomic operators Xa € A act as the iden- 
tity on the auxiliary space, i.e. / (g) Xa- Given these 
definitions, the derivation of the quantum filtering equa- 
tion remains essentially unchanged, so that the filter in 
either the operator form of ([s]) or the adjoint form of ^ 
is simply updated with the extended forms of operators 
given in the last paragraph. 

Since ^ is a classical parameter, we require that the 
reduced conditional density matrix {p^)t — Tr-^^ (pt) be 
diagonal in the basis of S. Thus we can write 



N 



(14) 



where 



Pi 



(0 - 



Tl■[m{^^\^I)Pt]^M\^^){L\®I] 

= E[i6)(c»i ® mm] = = mho.t])- (15) 



Then p^'-* is precisely the conditional probability for ^ to 

have the value S,i and the set {pj*' } gives the discrete con- 
ditional distribution of the random variable represented 
by S. Similarly, by requiring operators to be diagonal in 
Ti.^, we ensure that they correspond to classical random 
variables. In short, we have simply embedded filtering of 
a truly classical random variable in the quantum formal- 
ism. 

The fact that both states and operators are diagonal 
in the auxiliary space suggests using an ensemble form 
for filtering. As such, consider an ensemble consisting of 
a weighted set of N conditional atomic states, each state 
evolved under a different Later, in section IV we will 



call each ensemble member a quantum particle. For now, 
we explicitly write the conditional quantum state as 



N 



(16) 



1=1 



(i) . 

where is a density matrix on TIa- The reduced state, 
Tr-H^ (pf), is clearly diagonal in the basis of S. Using 
the extended version of the adjoint quantum filter in (|6| , 
one can derive the ensemble quantum filtering equations 



dpi 



i[^,Ho,p'i'>]dt+{Lp['^L^ - \L^Lpf^ - ^Pt 



dpf~> = (Tr[(L + Lt)pW] - Tr[/ ® {L + L^)pf] ) p\'^dW, 
dWt = dMt ~Tt[I^{L + L^)pf ] dt 



(17a) 



(17b) 
(17c) 
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(i) 

We see that each pi in the ensemble evolves under a 
quantum filter with H = S^iH^ and is coupled to other 
ensemble members through the innovation factor dWt, 
which depends on the ensemble expectation of the mea- 
surement observable. Note that one can incorporate any 
prior knowledge of ^ in the weights of the initial distri- 
bution {Po*'}. 

The reader should not be surprised that a similar ap- 
proach would work for estimating more than one param- 
eter at a time, such as three cartesian components of an 
applied magnetic field. One would introduce an auxil- 
iary space for each parameter and extend the operators 
in the obvious way. The ensemble filter would then be 
for a joint distribution over the multi-dimensional pa- 
rameter space. Similarly, one could use this formalism to 
distinguish initial states, rather than parameters which 
couple via the Hamiltonian. For example, in the case 
of state discrimination, one would introduce an auxiliary 
space which labels the possible input states, but does 
not play any role in the dynamics. The filtered weights 
would then be the probabilities to have been given a par- 
ticular initial state. In fact, using a slightly different 
derivation, Jacobs derived equations similar to (17) for 
the case of binary state discrimination |18| . Yanagisawa 
recently studied the general problem of retrodiction or 
"smoothing" of quantum states [25j . In light of his work 
and results in the following section, the retrodictive ca- 
pabilities of quantum filtering are very limited without 
significant prior knowledge or feedback. 



A. Conditions for convergence 

Although introducing the auxiliary parameter space 
does not change the derivation of the quantum filter, it 
is not clear how the initial uncertainty in the param- 
eter will impact the filter's ability to ultimately track 
to the correct value. Indeed, outside of anecdotal nu- 
merical evidence (which we will presently add to), there 
has been little formal consideration of the sensitivity of 
the quantum filter to the initial state estimate. Re- 
cently, van Handel presented a set of conditions which 
determine whether the quantum filter will asymptotically 
track to the correct state independently of the assumed 
initial state [12]. Since we have embedded parameter 
estimation in the state estimation framework, such sta- 
bility then determines whether the quantum filter can 
asymptotically track to the true parameter, i.e. whether 
linit Sij when ^ — ^i. In this section, we present 

van Handel's results in the context of our parameter es- 
timation formalism and present a simple check of asymp- 
totic convergence of the parameter estimate. We begin 
by reviewing the notions of absolute continuity and ob- 
servabihty. 

In the general stability problem, let pi be the true un- 
derlying state and p2 be the initial filter estimate. We 
say that pi is absolutely continuous with respect to p2, 
written pi <ti P2, if and only if ker pi D ker p2. In the con- 



text of parameter estimation, we assume that we know 
the initial atomic state exactly, so that pi <C P2 as long 
as the reduced states satisfy pf <C pf- Since these re- 
duced states are simply discrete probability distributions, 
and {{pI)2}, this is just the standard definition 
of absolute continuity in classical probability theory. In 
our case, the true state has (p|i^o)i ~ ^ij parameter 
has value ^i. Thus, as long as our estimate has non-zero 
weight on the i-th component, pi <ti P2- This is trivially 
satisfied if {p[''2q)2 for all j. 

The other condition for asymptotic convergence is that 
of observability. A system is observable if one can deter- 
mine the exact initial atomic state given the entire mea- 
surement record over the infinite time interval. Observ- 
ability is then akin to the ability to distinguish any pair 
of initial states on the basis of the measurement statistics 
alone. Recall the definition of the Lindblad generator in 
([3]) and further define the operator /C[X^] = L^^ Xa+XaL. 
Then according to Proposition 5.7 in [T^, the observable 
space O is defined as the smallest linear subspace of A 
containing the identity and which is invariant under the 
action of C and IC. The filter is observable if and only if 
A~0, or equivalently dim^ = dimC 

In the finite-dimensional case, van Handel presents an 
iterative procedure for constructing the observable space. 
Define the linear spaces Z„ C .4 as 



Zq — span{/} 

Zn = span{Z„_i,£[Z„_i],/C[2^„_i]}, n > 



(18) 



The procedure terminates when Z„ = Zn^^i, which is 
guaranteed for some finite n = m, as the dimension of 
Zn cannot exceed the dimension of the ambient space 
A. Moreover, the terminal Zm — O, so that using a 
Gram-Schmidt procedure, one can iterativcly find a ba- 
sis for O and easily compute its dimension. Note that 
for operators A and B, the inner-product {A, B) is the 
Hilbert-Schmidt inner product Tv[A^B]. 

Given these definitions, one has the following theorem 
for filter convergence and corollary for parameter estima- 
tion. 

Theorem 1. (Theorem 2.5 in fW^ ) Let Trf (X^) be the 
evolved filter estimate, initialized under state pi . If the 
system is observable and pi <C P2, the quantum filter is 
asymptotically stable in the sense that 



yXa e A 



(19) 



where the convergence is under the observations gener- 
ated by pi . 

One could use this theorem to directly check the sta- 
bility of the quantum filter for parameter estimation, us- 
ing the extended forms of operators in C and /C and 
being careful that the observability condition is now 
dimO — diml?(7ij) ® A. However, the following corol- 
lary relates the observability of the parameter filter to the 
observability of the related filter for a known parameter. 
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Combined with the discussion of extending the absolute 
continuity condition, this then gives a simple check for 
the stability of the parameter filter. 

Corollary 1.1. Consider a parameter ^ which takes on 
one of N distinct positive real values {^i}. // the quan- 
tum filter with known parameter is observable, then the 
corresponding extended filter for estimation of ^ is ob- 
servable. 

Proof. In order to satisfy the observability condition, we 
require dimO = Nr, where we have set dim^ = r and 
used the fact that dimS(7i^) = N. Given that the filter 
for a known parameter is observable, its observable space 
coincides with A and has an orthogonal operator basis 
{Ai}, where we take Aq = I. 

Similarly, consider the A'^-dimensional operator space 
D{TL^). If {^,;} are distinct, any set of the form 



Although these conditions provide a simple check, we 
would like to stress that they do not determine how 
quickly the convergence occurs, which will depend on the 
specifics of the problem at hand. Additionally, as posed, 
the question of observability is a binary one. One might 
expect that some unobservable systems are nonetheless 
"more observable" than others or simply that unobserv- 
able systems might still be useful for parameter estima- 
tion. Given the corollary above, one can see that this 
may occur if a single parameter =0. Then has a 
row of all zeros, so that the maximal dimension of a set 
of linearly independent powers of S is iV — 1 . Similarly, 
if one allows both positive and negative real-valued pa- 
rameters, the properties of are not as obvious, though 
in many circumstances, having both and — renders 
the system unobservable. We explore these nuances in 
numerical simulations presented in the following section. 



~.'"'},h£^,h^k.j\ii^ j (20) 



is linearly independent, since the corresponding general- 
ized Vandermonde matrix 



Si ?1 



(21) 



has linearly independent columns |26j . 

Following the iterative procedure, we construct the ob- 
servable space for the parameter estimation filter starting 
with /(g) Aq, which is the identity in the extended space. 
We then iteratively apply C and K, until we have an in- 
variant linear span of operators. The only non-trivial 
operator on the auxiliary space comes from the Hamil- 
tonian part of the Lindblad generator, which introduces 
higher and higher powers of the diagonal matrix S. Since 
dim£)(7i^) (g) .A is finite, this procedure must terminate. 
The resulting observable space can be decomposed into 
subspaces 



k] G N 



(22) 



where kl is some increasing sequence of non-negative in- 
tegers which correspond to the powers of S that are intro- 
duced via the Hamiltonian. Note that the specific values 
of kl depend on the commutator algebra of Hq and the 
atomic-space operator basis {Ai}. Regardless, since the 
Hamiltonian in C can always add more powers of S, the 
procedure will not terminate until Oi is composed of a 
largest linearly independent set of powers of S. This set 
has at most N distinct powers of S, since it cannot ex- 
ceed the dimension of the auxiliary space. Given that any 
collection of N powers of S is linearly independent, this 
means once we reach a set of N powers kl , the procedure 
terminates and dimO^ = N. Since O has r subspaces 
Oi, each of dimension N , dim O — Nr as desired and the 
observability condition is satisfied. □ 



Qubit Example 

Consider using the single qubit magnetometer of Sec- 
tion [n] as a probe for the magnetic field B. Since the 
initial state is restricted to the x-z plane, the y compo- 
nent of the Bloch vector is always zero and thus is not 
a relevant part of the atomic observable space, which is 
spanned by {J, Ca;, tiz}. In some sense, the filter with 
known B is trivially observable, since we assume the ini- 
tial state is known precisely. 

When B is unknown, the ensemble parameter filter is 
given by 



def^ = -2B,dt-t- KCOs(6»f^)(sin(6lf^) - 2{a^)^)dt 



+ 2-/kcosi 



s{ef^)dWt 



(23a) 
(23b) 



where dWt = dMt - 2VK(cr^)('^) and (ct^)-^ = 



'^iP^t^ sm{9f'). We simulated this filter by numerically 



integrating the quantum filter in ( 10 1 using a value for 



B uniformly chosen from the given ensemble of potential 
B values. This generates a measurement current dMt, 
which is then fed into the ensemble filter of ( 23 1 . For all 
simulations, we set k = 1 and used a simple Ito-Euler 
integrator with a step-size dt = 10~^ 

Figure |3ja) shows a simulation of a filter for the case 
B £ {2k, 5k, 8k, 12k}. The filter was initialized with a 

uniform distribution, p^ ^ = 1/4. For the particular tra- 
jectory shown, the true value of B was 2k and we see that 
the filter successfully tracks to the correct B value. This 
is not surprising, given that the potential values of B 
are positive and distinct, thus satisfying the convergence 
corollary. It is also interesting to note that the filter 
quickly discounts the probabilities for 8k, 12k, which are 
far from the true value. Conversely, the filter initially fa- 
vors the incorrect B = 5k value before honing in on the 
correct parameter value. 
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FIG. 3: (a) Filtered p^*' for B G {2k, 5k, 8k, 12k}. The filter tracks to the true underlying value of B = 2k (b) Filtered p'*' for 
B £ {— K, +k}. The filter does not track to i? = +k with probability one, though it is the most probable parameter value. 



In Figure p[b) , we show a simulation for the case of 
B e — k}, which does not satisfy the convergence 

corollary. In fact, using the iterative procedure, one finds 
the observable space is spanned by {I ^ I, I (g) (7^, B (g) 
ax,B'^ g) I,B^ (g) a:,,B^ g) a^}- But since B = {q \ ), 
B^ = K^I so that only 3 of the 6 operators are linearly 
independent. Although the filter does not converge to 
the true underlying value of -B = +k, it does reach a 
steady-state that weights the true value of B more heav- 
ily. Simulating 100 different trajectories for the filter, we 
observed 81 trials for which the final probabilities were 
weighted more heavily towards the true value of B. This 
confirms our intuition that the binary question of observ- 
ability does not entirely characterize the performance of 
the parameter filter. 

Figure |4] shows the rate of convergence of filters meant 
to distinguish different sets of B. The rate of convergence 
is defined as the ensemble average of the random variable 



1 , if pj'^ > a for any i 
0, otherwise 



(24) 



Although any individual run might fluctuate before con- 
verging to the underlying B value, the average of over 
many runs should give some sense of the rate at which 
these fluctuations die down. For the simulation shown, 
we set a = 0.95 and averaged /0.95 over 1000 runs for two 
different cases — either all possible B values are greater 
than K or all are less than k. As shown in the plot, the 
former case shows faster convergence since the B field 
drives the dynamics more strongly than the measurement 
process, which in turn makes the trajectories of different 
ensemble members more distinct. Of course, one cannot 
make the measurement strength too weak since we need 
to learn about the system evolution. 



IV. QUANTUM PARTICLE FILTER 

Abstractly, developing a parameter estimator in the 
continuous case is not very different than in the finite 
dimensional case. One can still introduce an auxiliary 
space Tij, which is now infinite dimensional. In this 
space, we embed the operator version of ^ as 



(25) 



where = $10 and {^\C) = <5(C - C)- Again, by ex- 
tending operators appropriately, the filters in ^ and ^ 
become optimal parameter estimation filters. We gener- 
alize the conditional ensemble state of ( 16 1 to 



(26) 



where pt{£,) = ^'(CI-^'^[o,t]) is the continuous conditional 
probability density. Although the quantum filter pro- 
vides an exact formula for the evolution of this density. 







«-B € {2k, 


4k, 6k} 
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FIG. 4: Rate of convergence (/0.95), averaged over 1000 tra- 
jectories. The filters are for cases when possible B values are 
either all larger or all smaller than the measurement strength 

K. 
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calculating it is impractical, as one cannot exactly rep- 
resent the continuous distribution on a computer. The 
obvious approximation is to discretize the space of pa- 
rameter values and then use the ensemble filter deter- 



mined by (17 1; indeed such an approach is very common 



in classical filtering theory and encompasses a broad set 
of Monte Carlo methods called particle filters piP, I^T] . 

The inspiration for particle filtering comes from noting 
that any distribution can be approximated by a weighted 
set of point masses or particles. In the quantum case, 
we introduce a quantum particle approximation of the 
conditional density in ( 26 1 as 



N 



MO 



(27) 



The approximation can be made arbitrarily accurate in 
the limit of ^ oo. Plugging this into (26), we recover 



precisely the form for the discrete conditional state given 
in ( 16 1 . Accordingly, the quantum particle filtering equa- 



tions are identical to those of the ensemble filter given in 
(17 1. The only distinction here is in the initial approxi- 



mation of the space of parameter values. Thus the basic 
quantum particle filter simply involves discretizing the 
parameter space, then integrating the filter according to 
the ensemble filtering equations. 

The basic particle filter suffers from a degeneracy prob- 
lem, in that all but a few particles may end up with neg- 
ligible weights pi . This problem is even more relevant 
when performing parameter estimation, since the set of 
possible values for ^ are fixed at the outset by the choice 
of discretization. Even if a region in parameter space 
has low weights, its particles take up computational re- 
sources, but contribute little to the estimate of More 
importantly, the ultimate precision of the parameter es- 
timate is inherently limited by the initial discretization; 
we can never have a particle whose parameter value S.i 
is any closer to the true value ^ than the closest initial 
discretized value. 

In order to circumvent these issues, we adopt the ker- 
nel resampling techniques of Liu and West [28j . The idea 
is to replace low weight particles with new ones concen- 
trated in high weight regions of parameter space. One 
first samples a source particle from the discrete distribu- 
tion given by the weights {Pt'^}, ensuring new particles 
come from more probable regions of parameter space. 
Given a source particle, we then create a child particle 
by sampling from a Gaussian kernel centered near the 
source particle. By repeating this procedure N times, we 
create a new set of particles which populate more proba- 
ble regions of parameter space. Over time, this adaptive 
procedure allows the filter to move away from unimpor- 
tant regions of parameter space and more finely explore 
the most probable parameter values. 

The details of the adaptive filter lie in parameterizing 
and sampling from the Gaussian kernel. Essentially, we 
are given a source particle, characterized by and 

(i) . 

p\ , and using the kernel, create a child particle, charac- 



terized by and pf' . One could attempt to sample 

from a multi-dimensional Gaussian over both the param- 
eter and atomic state components, but ensuring that the 
sampled p\ is a, valid atomic state would be non-trivial 
in general. There will be some cases, including the qubit 
example in the following section, where the atomic state 
is conveniently parameterized for Gaussian resampling. 
But for clarity in presenting the general filter, we will 
create a child particle with the same atomic state as the 
parent particle. 

Under this assumption, the Gaussian kernel for parent 
particle i is characterized by a mean /i^'^ and variance 
cr^^'\ both defined over the one dimensional parameter 
space. Rather than setting the mean of this kernel to 
the parameter value of the parent, Liu and West suggest 
setting 



M^'^ - <^ + (1 - o)l a e [0, 1] 



(28) 



where ^ — '^iP^t^S.i is the ensemble mean. The param- 
eter a is generally taken to be close to one and serves 
as a mean reverting factor. This is important because 
simply resampling from Gaussians centered at results 
in an overly dispersed ensemble relative to the parent 
ensemble. The kernel variance is set to 



2W 



h'Vu h e [0, 1] 



(29) 



where Vt — J2iPi^\^i is the ensemble variance and 

h is the smoothing parameter. It is generally a small 
number chosen to scale with TV, so as to control how 
much kernel sampling explores parameter space. While 
a and h can be chosen independently, Liu and West relate 
them by ft^ = 1 — a^, so that the new sample does not 
have an increased variance. 

Of course, it would be computationally inefficient to 
perform this resampling strategy at every timestep, es- 
pecially since there will be many steps where most par- 
ticles have non-negligible contributions to the parameter 
estimate. Instead, we should only resample if some un- 
desired level of degeneracy is reached. As discussed by 
Arulampalam et al. |21j . one measure of degeneracy is 
the effective sample size 



1 



(30) 



At each timestep, we then resample if the ratio Ncs/N 
is below some given threshold. We are not aware of an 
optimal threshold to chose in general, but the literature 
suggest 2/3 as a rule of thumb pU] . 

Altogether, the resampling quantum particle filter al- 
gorithm proceeds as follows: 



Initialization for i = 1, 



,N: 



1. Sample from the prior parameter distribu- 
tion. 
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2. Create a quantum particle with weight pf' = 
1/N, parameter state and atomic state 

Po "* — PO) where po is the known initial atomic 
state. 



Repeat for all time: 



1, 



Update the particle ensemble by integrating a 
timestep of the filter given in (17). 



If N^s/N is less than the target threshold, cre- 
ate a new particle ensemble: 

Resample for i — 1, . . . ,N: 

(a) Sample an index i from the discrete 
density {pt*"*}. 

(b) Sample a new parameter value from 
the Gaussian kernel with mean /i'*^ 




and variance (t^'*'' given by (|28| and 



(291 



FIG. 5: Kernel density reconstruction of pt{B)dB = 
P{B\M[o,t])dB for = 1000 particle filter set with dB = 
IOk/150, a = 0.98, h = 10"'^ and resampling threshold of 



(c) Add a quantum particle to the new 2/3. The true magnetic field was B = 5k. 



ensemble with weig ht pf^ = 1/N, pa- 
rameter state and atomic state 



Pt 



ii) 



(i) 
Pt 



Unfortunately, checking asymptotic convergence of the 
filter is more involved in the continuous-valued case, as 
the observability and absolute continuity conditions re- 
quire extra care in infinite dimensions. However, given 
that the quantum particle filter actually works on a dis- 
cretized space, in practice we can simply use the results 
we had for the finite-dimensional case. As before, we note 
that one can generalize the quantum particle filter to mul- 
tidimensional parameters by using a multi-dimensional 
Gaussian kernel. One might also consider using alternate 
kernel forms, such as a regular grid which has increasingly 
finer resolution with each resampling stage. We will not 
consider such extensions here. 



]ubit Example 



B range from to 10k in our units. The resampling 
parameters were a = 0.98, h = 10~^ and resampling 
threshold 2/3. Note that we chose not to use Liu and 
West's relation between a and h. 

In order to generate the figure, each particle's weight 
and parameter values were stored at 50 equally spaced 
times over the integration timespan. Using Matlab's 
ksdensity function, these samples were then used to re- 
construct Pt{B) via a Gaussian kernel density estimate 
of the distribution. The resulting kernel density esti- 
mate was then evaluated at 150 equally spaced B values 
in the range [0, 10k], which we plotted as pt{B)dB with 
dB — IOk/150. As is seen in the figure, after some ini- 
tial multi-modal distributions over parameter space, the 
filter hones in on the true value of B = 5k. For the sim- 
ulation shown, the final estimate was B — 5.03k with 
uncertainty ct^ = 0.18k. The filter resampled 7 times 
over the course of integration. 



We now consider a resampling quantum particle fil- 
ter for the qubit magnetometer introduced earlier in the 
paper. As hinted at in the previous section, since the 
qubit state is parameterized by the continuous variable 
Of, we can easily resample both the magnetic field Bi 
and state 9'^'^^ using a two-dimensional Gaussian kernel 
for (13,, 9'^'^), with mean vector and covariance matrix 



given by generalizations of ( 28 1 and ( 29 ) . Since different 



values of B result in different state evolutions, resampling 
both the state and magnetic field values should result in 
child particles that are closer to the true evolved state. 

Figure [5] shows a typical run of the quantum particle 
filter for N — 1000 particles. The true B value was 5k 
and the prior distribution over B was taken to be uniform 
over the interval [0, 10k]. As before, we used an Ito-Eulcr 
integrator with a step-size of dt = 10~^k. Note that both 
the timespan of integration and the potential values of 



V. CONCLUSION 

We have presented practical methods for single-shot 
parameter estimation via continuous quantum measure- 
ment. By embedding the parameter estimation problem 
in the standard quantum filtering problem, the optimal 
parameter filter is given by an extended form of the stan- 
dard quantum filtering equation. For parameters taking 
values in a finite set, we gave conditions for determining 
whether the parameter filter will asymptotically converge 
to the correct value. For parameters taking values from 
an infinite set, we introduced the quantum particle fil- 
ter as a computational tool for suboptimal estimation. 
Throughout, we presented numerical simulations of our 
methods using a single qubit magnetometer. 

Our techniques should generalize straightforwardly for 
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estimating time-dependent parameters and to a lesser 
extent, estimating initial state parameters. The binary 
state discrimination problem studied by Jacobs [TH] is 
one such example and his approach is essentially a spe- 
cial case of our ensemble parameter filter. We caution 
that the utility of initial state or parameter estimation 
depends heavily on the observability and absolute conti- 
nuity of the problem at hand. Future extensions of our 
work include exploring alternate resampling techniques 
for the quantum particle filter and developing feedback 
strategies for improving the parameter estimate. More 



broadly, we believe there is much to be learned from clas- 
sical control and parameter estimation theories. 
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